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Abstract 


In this paper we consider high-order centered finite difference approximations of hyperbolic 
conservation laws. We propose different ways of adding artificial viscosity to obtain sharp 
shock resolution. For the Riemann problem we give simple explicit formulas for obtaining 
stationary one- and two-point shocks. This can be done for any order of accuracy. It is 
shown that the addition of artificial viscosity is equivalent to ensuring the Lax £-shock 
condition. We also show numerical experiments that verify the theoretical results. 
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1 Introduction 


Centered finite difference methods, especially high-order ones, are often computationally 
efficient when the solution of the underlying problem is smooth. For non-smooth solu- 
tions, however, these methods produce excessive oscillations, which ultimately will ruin 
the solutions completely. One way to overcome the spurious phenomena is to introduce 
numerical viscosity which will smooth the numerical solution. But there are caveats; vis- 
cosity must not be used such that unnecessary smoothing occurs. As the computational 
mesh is refined the viscous effects should decrease while still damping the oscillations. 
Nessyahu and Tadmor [7] used a Lax-Friedrichs solver to construct non-oscillatory 2nd- 
order centered difference methods for hyperbolic conservation laws. A similar approach 
was undertaken by Harten and Lax in [3], where they constructed an approximate Rie- 
mann solver from a modified version of Richtmyer’s two-step scheme. In this paper we 
shall develop a general theory on how to achieve sharp shock resolution for high-order 
finite difference approximations of systems of conservation laws by adding artificial vis- 
cosity. We shall thus consider finite difference solutions of the Riemann problem 


u t + fx — 0 


u(x, 0) 


ul 

UR 


x < 0 
a; > 0 , 


(1) 


where it is assumed that the states ui , ur £ can be connected via a fc-shock moving 
with speed s\ f — f{u) £ 9^ is assumed to be differentiable. By means of the coordinate 

transformation 

y — x — st 
T — t 


the original Riemann problem (1) is transformed to a stationary problem 

(/ - ™)y = 0 

/ n x / u L y < 0 (2) 

u{y ’ 0) = \n R „> 0 , 

since u T = 0 in the (y, r)-coordinates. It will be assumed that the entropy solution of 
eq. (2) can be obtained as the pointwise limit (boundedly) of the regularized problem 

(f-su) v = L e u, lim u(y,<) = «L, lim u{y,t) = u R , (3) 

\ J y ii — <. — nr; V — ►CO 


when £ — > 0; L c is a linear elliptic operator. 


2 Second- Order Difference Methods 

When analyzing second-order difference methods we follow the principles set forth in [1 , 5]. 
Hence, eq. (3) is discretized as 

Do(fj — suj) = shD+D-Uj , £ > 0 , (4) 


1 



subject to the boundary conditions 

lim Uj = up , lim Uj — ur . (5) 

j — ►— oo j— +-oo 

The right-hand side of eq. (4) corresponds to choosing L e = edy , e — ► 0 + in eq. (3). The 
difference operators are defined as 

DoUj — 7^r ( w .?+ 1 j D~Uj w i“i) ’ 

where h = Xj+i — Xj is the uniform mesh size. The central difference Do(fj — suj) can be 
rewritten in conservative form as 

Do ( fj “ $ u j) = D+ ^-(/j — suj + fj-\ — ■ 

Thus, eq. (4) can be expressed as 

D + Q(/i ~ 5U j + /i-i ~ - e ( u i - u i-i)) = 0 » 

which in turn leads to 

2 (/j-J- 1 — su i + i + fj ~ su i) — e ( u j+i ~ u i) — 2 (fj ~ su i fj~ l ~ su j- 1) — £ ( u i ” u j-i) 


= • • • = fL - SU L , 


where the last equality follows from the boundary conditions (5). Consequently, 
fj+ 1 - /l - a(u,-+i - «l) - 2eu J+1 = -(/,• - fp) + s(uj - up) - 2 eu 3 . 
Let 

F(u) = f(u) - fp - s(u - up) -2 eu 
G(u) = f L - f(u) + s(u - u L ) -2eu. 

The difference approximation (4), (5) can then be written as 


( 6 ) 


( 7 ) 


Letting j —* oo yields 


F(«r) = G(u h ), 


i. e., 


/fi — /l = — • 


( 8 ) 


This is the familiar Rankine-Hugoniot condition, which is fulfilled since up and ur are 
the states on either side of the stationary fc-shock problem (2). We can thus define F(u) 
and G(u) as 

F(u) = /(«) - fp - s(u - up) - 2 eu , , 

G(u) = fp — f(u) -{- s(u — up) — 2su , 

where up = up or up = ur. The second choice corresponds to letting j — *• oo in eq. (6). 
The Rankine-Hugoniot condition would then follow by taking the limit j — * —oo in eq. (7). 
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2.1 Nonlinear Analysis of One- and Two-Point Shocks 

Given that there exists a solution, eq. (7) can be interpreted in the following way: The 
two states u R and u R can be connected to each other via a stationary viscous profile 
{ujlToo where Uj satisfies the nonlinear recursion formula (7). The question then arises 
naturally whether it is possible to connect two given states u L and u R by a viscous profile 
consisting of a finite number of intermediate states um,, j = 1 !•••>?• We begin by proving 
a negative result. 

Proposition 2.1 There are no states u L / u R such that u L and u R can be coupled by 
the viscous regularization (7) without any intermediate state , i. e., there are no zero-point 
shocks. 


Proof: Suppose that the two states u L and u R can be directly coupled by eq. (7), l. e., 

F(u r ) = G(u l ). 

From definition (9) it follows immediately that 

Jr — fp ~ s{ u R — u p ) — 2eu.R = fp — /z, + ${ul — u p) ~ 2 eu R . 

But the Rankine-Hugoniot condition (8) implies that 

2 eu R = 2 sul . 

Since £ > 0 it follows that u R = u L , which proves the proposition. D 

Proposition 2.2 Two states u L , u R G $ d can be connected via a viscous profile (7) using 

one intermediate state u M € (one-point shock). Furthermore, if the states u L and u R 
are close, the Lax k-shock condition is equivalent to having £ > 0 in eq. (4)- 

Proof: A single intermediate state u M must according to eq. (7) satisfy 

G(u l ) = F(u m ) (io) 

G(um) = F(u r ) . 

The definition of F and G and the Rankine-Hugoniot condition (8) together imply that 

-2eu L = Zm - Jl - s(um - U L ) - 2eu M m) 

-2 eu R = Jl - Im + ~ u l ) - 2 £um 

Adding these two equation yields 

UM = + u fl) • 
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We have thus found the intermediate state um- The scalar viscosity e and the shock 
speed s remain to be determined. Using the well-known Roe-linearization one can re- 
write eq. (11) a-s 


um) i^L u m) — 2 E{uij um) 

( A{um , Ur) — si) (um — Ur) = —2 s(um ~ Ur) , 
where the Roe-matrices A(ul,um) and A(um,ur) are defined by 

A(u l ,u m )=[ f'(0u L + (\ -d)u M )dO 

Jo 

a(u m ,ur)= f f'(eu M + (\-e)u R )d0. 

Jo 

But eq. (13) has a solution if and only if 

2e = Hk(u L ,u M ) - s , 2e = -Hk{uM,UR) + s , 


(13) 


i. e., 


£ = -{Vk{u L ,U M ) ~ flk(uM,UR )) 


(14) 


•s = -(fik(uL, Um) + Pk(uM, Ur)) . 

Here hic(ul,um) and Hk{uM,up) denote eigenvalues of A(ul,um) and A(umiUr), respec- 
tively. From the definitions of the Roe-matrices it follows immediately that 

A(u L ,U L ) = f'{u L ) , A(uh, Ur) = /'(ur) . 

Denote by \k(u) the eigenvalues of f'(u). Hence, Taylor expansion yields 

Pk{uL,UM) = Ajt(ut) + 0(\um - Ul|) = A* (ui) + 0(\ur - Ul\) , 

where the second equality follows from eq. (12). Similarly, 

H( u M, Ur) = A k (u R ) -)- 0(\ur - U L |) . 

The Lax fc-shock condition requires that 

A k{uL) > S > A k(u R ) 

hold for exactly one characteristic family k, \ < k < d. Thus, for weak shocks the Lax 
fc-shock condition holds iff 


Vk{uL, Um) > S > Hk(uM , Ur) , 
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which in turn implies that £ > 0, i. e., eq. (4) satisfies an entropy inequality. Conversely, 
suppose that e > 0. Then 

Hk{uL,u M ) > 


Hence, 

Hk(uL, u m) > -(^k{ u Li u m) + fJ-k{UM, u R)) > ^k{uMi u R) ■ 

Recalling eq. (14) we see that 

Hk(uL,UM) > s > Hk{uM,UR). 

Thus, for weak shocks our solutions obey the Lax fc-shock condition. 

Finally, we need to verify that s defined by eq. (14) is compatible with the Rankine- 
Hugoniot condition (8). We have 

f L - f R = A(u l , u m ){u l - u M ) + A(um, ur)(u m ~ ur) . 

For a one-point shock, however, ul — um, um — U R are eigenvectors of A{ul,um) and 
A(U M ,ur) (cf. eq. (13)), respectively, whence 

fl — Ir = lik{uLi Um){ul — Um) + ^k{ u Mi Ur)(um — Ur) . 


Using eq. (12) yields 

h - }r - U m) + Hk{uMi Wfi))(“L - ur) . 

Thus, the Rankine-Hugoniot condition (8) is satisfied with 


s = -( Hk{uL , u m) + ur)) . 

This concludes the proof. ^ 

Remark: The proof of the previous proposition shows that one should choose e according 
to eq. (14). In the actual implementation, however, it would be advantageous to use 

£ = -rl^kiuL, u M ) ~ f*k(uM,UR)\ , 

since the argument of the modulus function would be positive for a true fc-shock. If, on 
the other hand, the argument should happen to be negative due to round-off errors, for 
instance, then the modulus function will prevent the formation of an entropy violating 
shock. ^ 


In the scalar case £ and s can be expressed directly in terms of f(u) since 

h ~ Im 2 (/l - Im ) 


Hk{u L ,u M ) = A(u l ,u m ) = 


Ut — Um 


ut — u R 


fik { uM , u R ) = A(um,ur) = 


l m — Jr 

um — U R 


2(/m ~ /fl) 

UL — UR 
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Using these expressions in eq. (14) yields 

h — 2 Jm_ + Ir 

2 (u L - Ufl) 

II ~ Ir 


£ — 


( 15 ) 


3 = 


UL - U R 


The next proposition shows that the smallness assumption of proposition 2.2 is not needed 
for strictly convex scalar functions f{u). 

Proposition 2.3 Let f(u) be a twice differentiable , strictly convex scalar function, and 
let £ be given by eq. (15). Then 

e > 0 f{u L ) > f'(u R ) . 


Remark: The latter condition is the well-known entropy condition for scalar conservation 
laws [6]. □ 

Proof: Taylor expansion yields 

Jr = f(u R - u M + U M ) = f{u M ) + f(u M ){uR - U M ) + i f"(w)(u R - u M ) 2 

fh = f(u L -u M + u M ) = f(u M ) + f(u M )(uL ~ U M ) + t/"( v )(u L - u M f 
for some v and w . Using eq. (12) then gives 

fh - 2/m + Ir 1 , ,// 


£ = 


2(ul — Ur) =16 


But /" > 0 because of strict convexity. Hence, 

£ > 0 •<=> U L > UR «=> f(u L ) > /'( Ur ) . 

The second equivalence follows since f(u) is a strictly increasing function of u (because 
of strict convexity). □ 

Example: For Burgers’ equation the flux is given by 

1 


From eq. (15) it then follows that 


/(«) = ^ . 


e - 1(«L - Ur) 

•s - 1(«L + Ur) . 


(16) 
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□ 


In case of a convex scalar function f(u) the one-point shock condition (10) can be 
depicted as (recall that F(ul) = G{ui ), F(ur) = G{ur)) 



ourselves to one-point shocks, which the following discussion will show. Suppose that we 
want to connect ul and u R by a viscous profile using two intermediate states um 2 and 
um 2 - This amounts to requiring 

G(u l ) = F(u Ml ) 

G(u M i ) = F(um 2 ) 

G{um 2 ) = F(ur) . 

Adding these three equations yields 

um 2 - «Mi = \( UR ~ UL ) ‘ 

As before we rewrite the first and third equations of (17) as 

(A(ul, UMi) - si) ( UL - UMi) = 2 e(ul — UMj) 

(A(um 2 ,Ur) - si) ( Um 2 - Ur) = - 2 e(um 2 - Ur) , 

which can be solved iff 

2e = (i k {uLiU Ml ) - s > 0 u Ml ur — a^r k (u L ,u Ml ) 

-2e = Vk(uM 2 , ur) - s < 0 u M2 - ur = a 2 r k (uM 2 ,UR ) , 

where r fc (u L ,u Ml ), r k {u M2 ,u R ) denote the fcth eigenvectors of A{u L ,u Ml ), A(u M2 ,u R ) 
The inequalities are immediate consequences of the A:-shock condition. The first set of 


(17) 

(18) 
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equations yields 

e = ~(n k (u L ,u Ml ) - Hk{uM 2 ,u R )) 

( 19 ) 

s = 2^ k ^ UL,UM ^ f l k{ u M 2 ,U R )) • 

In particular, if and um 2 are known, so are e and s. The intermediate states are 
determined by eq. (18), by 


u Mi U L — ^-X^kiuLt ) (201 

um 2 - Ur - otir k (uM 2 ,u R ) , 1 

and by the Rankine-Hugoniot conditions (8). As opposed to the case of a one-point shock, 
the explicit eigenvector structure of A(ul,um^) and A(um 2 ,u r ) must be known in order 
to compute the intermediate states. 

Example: In the scalar case eq. (20) becomes redundant (rjt = 1). We are thus left with 


um 2 ~ wjw, = {u R — ul)/2 

/r — ft = s(u R — Ul ) 

s = {fi(u L ,u Ml ) + h(um 2 ,u r ))/ 2, 


[l - fau + hh - [r _ 2(/l - /r) 
ul ~ umi um 2 — u R ul — u R 

um 2 ~ u Ml = ( u R - u L )/ 2 . 

For Burgers’ equation one gets the linear system 

umi + u m 2 = u l + u R 

UMi — um 2 = ( ul — u r )/2 , 

which yields 

u Ml = (3u l + u R )j 4 
um 2 = (ul + 3uh)/4 . 

Consequently, 

M( u L, Um,) = g(? U L + U r) 1 t J '( u M 2 ,U R ) = -(ul + 7 U R ) , 

which implies 

£ = - «r) 

5 = l( u z. + u fi) • 


( 21 ) 


( 22 ) 


□ 
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2.2 Linear Perturbation Analysis 

In general we cannot expect the methods based on nonlinear analysis to be insensitive to 
perturbations. A different approach could be based upon linearization of eq. (7) [1,5]. 

{f (up) — sl)uj . |-i — 2euj+i = — (/ (up) — sl)uj — 2 euj . 

This equation can be diagonalized when f'(u) corresponds to a hyperbolic operator. The 
diagonalization is formally obtained by letting f'(u P ) -* \p, I -+ 1. We thus obtain 

(A p — s — 2e)u J+] = ( A p T s 2 c)uj . 

Suppose that u P = u R . This corresponds to linearizing around the state to the right of 
the shock. We can thus express uj+ 1 as a function of uy 

— A fi -f 5 — 2s 
Uj+1 = Ah - 7 - 2 £ U > ’ 

where Uj is assumed to be to the right of the shock. Note that one could not reverse 
the recursion above, since one would ultimately cross the discontinuity, across which the 
linearization has no meaning. The linearization implies that Uj is viewed as a perturbation 
around the constant state ur. No matter what the value of Uj is, we can make the 
perturbation disappear in the next step by setting 

s — 2e — \r. ( 23 ) 


Similarly, linearizing around the left state up gives 

\p — s — 2e 

J -X L + s-2e 

where it now is assumed that Uj+i is to the left of the shock. Again, requiring 


s -f - 2e — \p 


(24) 


implies that perturbations are annihilated in one step. Combining eqs. (23), (24) yields 


e = -(A l - ^h) 

(25) 

5 = ^(A l + Ah). 

The above expression for s states that s is the arithmetic average of the characteristic 
speeds on either side of the shock. This is the correct value modulo second order terms 
[6]. In particular, one should expect the corresponding numerical scheme to work well 
for weak shocks. Note that the analysis thus far - linear as well as nonlinear - has been 
based on eq. (7), which was obtained from eq. (4) by factoring out the operator D+. 
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Example: For Burgers’ equation we have A p = up, whence 

£ = ~ u R ) 

2 (26) 
•s = ~{u L + Ur). 

□ 

In the above example the linear approach resulted in twice as much viscosity as the one- 
point-shock. Note that the shock speed s is unchanged. This is, however, a coincidence. 
The values of e and s are determined by two criteria, namely that perturbations to the left 
and to the right of the shock are annihilated in one step. In general, these requirements 
are incompatible with correct shock speed (cf. eqs. (14), (25)), i. e., the Rankine-Hugoniot 
condition. Of course, if we confine ourselves to eliminating the oscillations on only one 
side of the shock, then we can use the correct shock speed. 

The transformation of eq. (1) to the time independent problem (2) was done to enable 
the theoretical analysis. In practice, one often computes in the fixed coordinate frame of 
eq. (1), which requires no a priori knowledge about s; correct shock speed will follow from 
the conservation form of (1). However, one can use the values for e obtained from the 
theoretical analysis and still obtain good results, cf. subsequent sections. Furthermore, in 
practical implementations it would be desirable to implement the viscosity locally around 
the shock. This can be done by introducing a switch so as to turn off the viscosity in the 
smooth regime. One way to do this would be to replace the right-hand side of eq. (4) by 

eD+r J _i/ 2 D-Uj , 

where rj = 1 close to the shock, rj = 0 otherwise; is the interpolation of r 3 at the 

cell interface £j_i/ 2 ; ui and Ur would be replaced by some interpolated value of u to the 
left and the right of the shock, respectively. 


3 High- Order Difference Methods 

We shall now use the results from the previous sections for 2nd-order methods to construct 
artificial viscosity for high-order methods. The standard explicit centered approximation 
of djdx of order 2 r has the form 

Qlr — RirDo, (27) 

R-2t = £(-l Ta v {h 2 D+D-Y , (28) 

i /=0 

where the coefficients a„ are defined by 

c*o = 1 , 

^ = 4772^-1 » v = l, 2 ,...,r- 1 . 


(29) 
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We can also use compact implicit difference approximations of Pade type. Some of 
these operators can be written in the form (27), but with a different operator R 2r . For 
example, the compact 4th-order approximation is (27) with 

R 4 = F$ ) = (I + jD + D.)~ 1 . ( 30 ) 

We will always assume that the implicit operators have the form /?F) P , where P is 
a nonsingular explicit operator of finite bandwidth. 

3.1 Factorization of High-Order Dissipation 

The idea to use to construct the viscosity, is to use high-order dissipation heD , where D 
can be factored as D = R 2r D + D The approximation for the transformed equation (2) 

IS 

R 2r D 0 {f - av)j = heR 2r D + D-Uj , ( 31 ) 

where £ is a parameter. The boundary conditions are as before 

lim Uj = u L , lim Uj = u R . (32) 

j— ►— oo 3 * co 

By (31) we see that it is natural to consider the operator R 2t in the space M of grid 
functions {u,} with 

lim Vj — 0 . (33) 

j-*±oo J 

By definition, the inverse of the implicit operators exist. For the explicit operators we 
have 

Lemma 3.1 The explicit operator R 2r in (28) is non-singular in the space M. 

Proof: Consider the equation 

R-lrVj = 0 , j = 0 , ± 1 , • • • ( 34 ) 

for real vj, and let the scalar product and norm be defined by 

OO 

(t7,tl>)= E VjWjh, IMI 2 = (v,v) . 

j = - OO 

Summation by parts yields (using (34) and (33)) 

0 = (v, R 2r v) = ||w || 2 + (■ h 2 D+D.yv ) 

V=1 


= llvll* + S . 


1/=1 
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Note that the boundary terms disappear because of (33). Since all a u are positive, it 
follows that v — 0, which proves the lemma. □ 

Since R 2r is non-singular, the equation (31) holds if and only if 

D 0 (f — su)j = 1ieD+D-u 3 (35) 

holds. But this is the three-point scheme that we have analyzed in section 2. Therefore, 
the high-order methods with artificial viscosity, as in (31), produce exactly the same solu- 
tion as the three-point scheme does. To take advantage of the high-order approximation 
property, we must obviously implement the viscosity locally around the discontinuity. We 
use the same switch function r 3 as for the three-point scheme. We summarize the results 
in 

Proposition 3.1 The propositions in section 2 concerning zero- and one-point shocks for 
2 nd-order methods apply to high-order methods of the form (31) as well For one-point 
shocks e should be chosen according to eq. (14)> Here R 2r is defined either by the explicit 
formula (28 ) , or by an implicit operator as described above. □ 


In section 2 we also determined the viscosity coefficient e based on linear analysis. 
The linearized equation for the high-order methods is 


R 2 rDo(f , (up) — sl)uj = heR 2T D+D-Uj . (36) 

The arguments in section 2 can be applied to this equation, and the optimal choice of e 
for the three-point scheme is optimal in the same sense for the high-order method. 

For explicit operators R 2t , the effect of the artificial viscosity can be interpreted in 
the following way. Upon diagonalization, the general solution of (36) is given by 

«i = £ , (37) 

J/=l 

where k v = ^(e) are the roots of the characteristic equation; a u are arbitrary scalar 
coefficients. For e = 0 there are two roots 

Ki = 1 , /v 2 = — 1 - 

There are also r — 1 roots {ac^}^ 1 with (^(0)1 < 1, and r— 1 roots {k u }^ 2 with | «*,(()) | > 1, 
see [9]. All roots except K\ — 1 give rise to parasitic solutions, and k 2 — — 1 is the one 
that causes the trouble. The remaining roots also induce errors, but they are less severe. 
If \k u \ > 1, then its presence near the shock is not felt, since the solution is bounded 
as j — ► oo. The analogous arguments hold for the solution to the left of the shock. 
The special choice (25) of e for the three-point scheme gives k 2 = 0. The factored 
form (36) implies that the coefficients ov, v — 3,4, ...,2r vanish. This follows since 
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/?2r is nonsingular on the space Af, which means that eq. (36) is equivalent to eq. (4). 
Consequently, the general solution (37) is identical to that of a three-point scheme. But 
that is possible iff the coefficients <j v — 0, v = 3, 4, , 2 r; the roots v = 3, 4, . . . , 2r 
are non-zero, but they do not influence the solution. Hence, the solution of the linear 
equation has the form 

u, = <7, , (38) 

where <j\ is determined by the condition at xj = oc, i. e., ct\ = u#. 


3.2 High-Order Dissipation Based on Perturbation Analysis 

We shall now consider a different form of artificial viscosity for the explicit approximation. 
For the nonlinear equation we use 

RirDotf - su)j = j2 h^-'e^D+D-Yuj . (39) 

l/=l 

This choice of viscosity corresponds to having L e = YX=~i £ ^y U i £ » ~ * ® * n e 9- (3)- There 
are r viscosity coefficients to determine, but the width of the total difference operator is 
not increased by the viscosity terms. The linearized equation is 

R 2 rDo(f’(u P ) - sl)u j = jz h 2v -'e„{D+D-Yu 3 . ( 40 ) 

v = 1 

The general form of the solution is still given by (37). Instead of forcing the coefficients 
cr 3 , (T4 , . . . , cr 2r to vanish, we now choose e„ such that 

K 2 = ^3 = ... K r +1 : 0 , 

where, as before, 

|/e„| < 1 , v — 3, 4, ... r + 1 

when no viscosity is present. This procedure can also be viewed as reducing the linear 
approximation to an r -f- 1 point scheme near the shock. As above we implement the 
approximation with a switch r 3 . After a re-normalization (j -*■ j - N) the solution of the 
linearized equation can be expressed as 

Uj = <Ti + ’ l^l > 1 • 

u—r-{-2 

The parasitic part of the solution represented by the sum does not cause any harm. Even 
if the computation is carried out over a finite domain, such that j < N, where N is 
fixed, the coefficients a u are bounded since the solution is bounded at j = N . Therefore, 
< 7 I ,k{~ n is small near the shock where j « N. We do not expect any difficulties, even if 
the shock passes through the boundary, since our methods are either strictly or strongly 
stable, see [8, 2]. 
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In the analysis we have assumed that all roots k u are simple. Should there be multiple 
roots, can always be partitioned into one group with exactly r — 1 members inside 

the unit circle and another group with r — 1 members outside the unit circle. Hence, our 
principle still applies. All the tc u inside the unit circle vanish by our choice of e v . 

We illustrate our principle by considering a 4th-order approximation Q 4 . After diag- 
onalization eq. (40) becomes 


(12£2 + A p — s)uj+2 "I" (12ei — 48£2 — 8(Ap — s))uj+j — (24 £j — 72£2)u^ 
+(12 £i — 48£2 T 8(Ap — s))uj_i + (12£2 — A p + s)tij_2 = 0 


(42) 


Suppose that we have linearized eq. (39) to the right of the shock, i. e., up = ur. We now 
want to eliminate the characteristic roots «2 = — 1 and /C3, |/t3 1 < 1. Being to the right 
of the shock this corresponds to setting the coefficients in front of Uj_ 1 and ttj_2 to zero 
(j = A:, k + 1 , . . .). Hence, 

£1 = -^(A R - s) 

(43) 

e 2 = A(A R - s ). 

On the other hand, linearizing eq. (39) to the left of the shock corresponds to setting 
up = up in eq. (42). This time we set the coefficients in front of Uj+i and ttj+ 2 to zero 
(j = — fc, —k — 1 , . . .), thus resulting in 


£l = g(A L-S) 

= -^(A L-s). 


(44) 


One realizes immediately from eqs. (43) and (44) that £1 > 0, £2 < 0 is equivalent to the 
A;-shock condition A 1 > s > \r. The viscosity coefficients will be uniquely determined iff 
(cf. remarks at the end of section 2.2) 

•s = • 


Thus, 


£ 1 = t(Al - a r) 
o 


£ ! = _ •'») 


( 45 ) 


s — 


^(Ai + Afi) ■ 
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( 46 ) 


Using these values in eq. (42) yields 

(A^ — Xp)u J+2 — 8(Al — \p)u }+ i + 7(A l — A R)uj 

-8(A p - Xr)uj-i + (A F - X R )uj- 2 = 0 . 

Hence, to the right of the shock we obtain (A L ± X R ) 

■Uj+2 — 8«j+i + luj = 0, j = k, k + 1, . . . 

Similarly, 

7 uj — 8uj_i 4- tij - 2 — 0 , j — k , k 1, . . . 

to the left of the shock. Using the substitution j — ► —j in the latter equation and defining 
vj = u_j, the two recursive expressions above will coalesce into 

Uj + 2 — 8u ; +i + 7uj = 0 , j = k, k + 1, . . . 

whose characteristic roots are given by 

Kj = 1 , «4 = 7 . 

Summing up, choosing s, £i, and e 2 according to eq. (45) in eq. (40) (for r = 2) yields 
the following characteristic roots 

Ki = 1 , k 2 = 0 , «3 = 0 , « 4 = 7 . 

Thus, the parasitic modes corresponding to k 2 and « 3 have been eliminated. 

Example: For Burgers’ equation we have A p = up, whence 

£i = c{ul - ur) 
b 

e 2 = ~^ UL ~ UR ) 

5 = “(UL + Ur) . 

□ 


Numerical results for this approximation will be presented in section 5. 


3.3 Averaging Operators 


Yet a different kind of viscosity, based on simple 2nd-order averaging, was used in [2]. We 
consider the approximation of the problem in its original time dependent form 




(48) 
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Let 


s" +1 = p( fi ") 

be a time discretization of Runge-Kutta type of the high-order semi-discrete difference 
scheme (48). Then the algorithm is 


6“ +I = PK) 

«” +1 = (i"-l + 2i” + ' + aji')/4 

Linearize around a constant state away from the shock as was done earlier. Let 


(49) 


-4^ = S 2r Vj , S 2r = -Q 2T f(up) 

be the linearized difference approximation in space without any viscosity. Then the Runge- 
Kutta method is 

v? +l = P(kS 2r )v] , (50) 

where 

m 

P(z) = l + E^ 

is=l 

is a polynomial (third degree in our experiments), and k is the time step. The averaging 
procedure can be written as 

v” +1 = (/+ jD + £>-)5" +1 . 

Hence, one complete step with the filtered Runge-Kutta method for the linearized equation 
is 

vj +1 = P(kS 2T )Vj + jD + D.P(kS 2T )v 

For explicit approximations Q^ r the artificial viscosity contains difference operators of all 
orders between 2 and m(2r — 1) -f 2, i. e., odd order terms will appear. For implicit 
approximations Q^r all grid points are involved when the filter is applied at any given 
point. With steady state solutions one obtains 

\{I - P(kS 2r )) Vj = ~D + D.P(kS 2r ) Vj , 


where A = k/h. Hence, the viscosity coefficients depend on the Courant number. In the 
last section it is demonstrated that the averaging technique (implemented with a switch 
function) is also very efficient for non-stationary shocks. 
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4 Linear Discontinuities 


( 51 ) 


Consider the linear problem 

Uf fAUx H 

u(z, 0 ) = <~p(x) 

where y is the (possibly discontinuous) periodic initial data; // € 3R. The equation above 
can be discretized as 

(52) 


u 


n+1 = U n } - kfiDou ] + ekhD+D.u ” , 


u" +1 = A ( £ + f) U ? -1 + ^ ~ A ( £ ~ 2) 

where A = k/h. It is well-known that by choosing 

\jA 

e 2 

one gets an upwind scheme. Furthermore, requiring 


u 


j + 1 1 


A 2e |/x| 


*_± 

M 


(53) 


(54) 


results in the method of characteristics, i. e., 

u” +1 = if // > 0, u" +1 = u" +1 if fi < 0 . 

Now suppose that d x is approximated using the high-order operator Q 2r = R 2r D Q . We 
can then use the idea from section 3.1 to obtain the following implicit scheme 


f? 2r u" +1 = R 2r u™ ~ kfj,Q 2r u'j + ekhR 2r D + D-u r - . (55) 

Assuming that there exists a periodic solution of eq. (55), it follows that R 2r exists. In the 
class of periodic solutions eq. (55) is thus equivalent to eq. (52). Hence, the conclusions 
above also pertain to eq. (55). 

In regions where the flow is smooth it should not be necessary to resort to the method 
of characteristics in order to compute solutions of high accuracy. Consider an explicit 
4th-order method for convenience. The viscosity term is then modified according to 

h 2 

shR 4 D+r 3 _ 4/2 D- , R 4 = I - -D+r 3 _ 1/2 D- , (56) 

where r 3 _ 1/2 = 0.5(r j _i + r 3 ) is the interpolant of a switch function r 3 ; Q 4 = (I - 
(/i 2 /6) D+D-) D 0 - We have used a switch proposed by Jameson [4] 

( |A+u ; ~ A_Uj| V (57) 

~ \|A+Uj| + |A_Uj|J ' 
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Note that r 3 = 1 if A and A+u, do not have the same sign. In particular, r, = 1 
where there are high-frequent oscillations. Choosing p — oo ensures that r, = 0 at the grid 
points where the flow is smooth. Thus, in smooth regions we have R 4 = I. Assuming that 
the spurious phenomena are localized to a compact set it follows that R 4 is non-singular. 
Eq. (51) is now discretized by 


R 4 u^ = i? 4 u" — A T 4 u" 

R 4 uf^ = aj/? 4 u” + fijRiU'p — 7 j AL 4 u* 1) 
# 4 u" +1 = 6ji? 4 n" + p J R 4 u'P — QXL 4 u , 

where # 4 is defined by (56); the spatial operator L 4 is given by 

L 4 — hQ 4 — sh? R 4 D -|_r j—i /^D _ 


and the integration parameters are given by 

Oj = 3/4 

h = i/4 

7j = 1/4 
6, = 1/3 
V: = 2/3 
0 = 2/3 


if r j = 0 , 


aj = 0 
& = 1 
7i = 0 

Aj = 0 

Vi = 1 
Cj = 0 


if rj = 1 . 


(58) 


(59) 


(60) 


Now, in smooth regions where r 3 == 0 the discretization (58) reduces to (i? 4 = /, L 4 = Q 4 ) 

uj 1 * — Uj- \Q 4 u " 



+ 1’4 ,) 



( n+l _ 


= -u" -| ti, 

<1 J ' Q J 


( 2 ) 


2A ^ (2) 


which is a 3rd-order TVD Runge-Kutta method. Note in particular that the value of a is 
irrelevant whenever r 3 — 0. When r 3 = 1, on the other hand, we recover 


— i? 4 tij — fc/iQ 4 t£j -|- D—U™ . 


Should r j — 1 for all j this is nothing but the Euler forward method (where we have 
assumed periodicity) that reduces to an upwind method for a = |/z | /2 and to the method 
of characteristics for A = l/|/x|. If r 3 — 1 locally, we expect (58) to behave roughly like 
an upwind scheme in that region. 


Returning to the original formulation (58) we observe that since J? 4 is non-singular, 
eq. (58) can be written in explicit form 

uf = Uj — Xv n } 

uf* = otju" + PjU ( p - 7 J Avj 1) (61) 

u" +1 = SjU 1 - + rjj u j 2) - , 
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where v n , t> (1) and v (2) are the solutions of the tridiagonal systems 

R 4 v] = L 4 u] , R 4 vf = L A uf , * = 1,2. 

The tridiagonal structure of R 4 is given by 

R 4 Vj = \ (-rj-i/ 2 Vj-i + (6 + r J+1 / 2 + rj_i/ 2 )uj - r J+ i/ 2 Uj + i) • 

D v 

In the next section we will present numerical results obtained from (61). 


5 Numerical Results 

We begin by studying the factorized artificial viscosity described in section 3.1. We have 
used the explicit and implicit 4th-order schemes (31), obtained by setting 

r 2 , = R, = / - jD + D . , 4 ? = 4 ° = (/ + j D + D -~r ' . 

for solving Burgers’ equation (f(u) = u 2 / 2) with a stationary shock (s = 0). An explicit 
3rd-order Runge-Kutta method has been used to solve a time dependent problem in a fixed 
coordinate system. Thus, the shock speed s does not appear explicitly in the equations. 
As initial data we have taken u(x,0) = -x, which should result in a stationary shock at 
the origin for t > 1. Furthermore, u L = 1 and = -1. From proposition 3.1 it follows 
that choosing e according to eq. (16), i. e., e = 1/4, should result in a one-point shock. 


£ 

+• 

i 

I 


Fig. 1: One-Point Shock, 4th-order Explicit (+) and Implicit (o), £ - 1/4 

Below is the result when the viscous terms are turned on only in a neighborhood of 
the shock. The viscosity is turned on and off by the switch defined by eq. (57). We also 
verify the theoretical results for the two-point shock. The coefficient E is then given by 
eq. (22), that is, e — 3/8. 


u(x,0)--x. u{- 1 ,t) - 1 . u(1,t)--1, t-2, n-101 
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oz 


X}iso3stj\ peooq ‘g/g = 3 ‘(o) ^piidmj puB (-[-) ipijdxg aapjo-ip^ ‘qooqg ^uioj-om^ 'p '3ig 



001 -U 2-1 •I'-ft'On l * (l‘l*)n 'x--{o‘x)n 

8/C = 3 ‘(°) ^PH dui I P UB ( + ) !P!I dx a japjo-ip^ ‘qooqg iuioj-omx :g '2y 



00 t-u ' 2-1 ’i*“(n)n ‘l-(n-)n ’x--( 0 ’x)n 

X^isoosiyY i^ooq f/\ = 9 { (o) iioiiduij pire (-f) ipijdxg ispjo-ipp ( )poqg : Z Sij 







We now proceed to case where the shock is non-stationary. Again we compute in a 
fixed coordinate system using the factorized viscosity of section 3.1. Thus, the theory 
developed in the previous sections does not apply directly. It is still interesting to see how 
well this principle of introducing viscosity fares in practice. Indeed, it turns out that the 
viscosity coefficient s obtained through linearization (26) works better in this case. In all 
of the following numerical experiments we have used the adaptive switch (57) to turn on 
viscosity locally around the shock; u(x,0) = x < 0, and u(x,0) = ur, x > 0; ui = 2, 
u R — 0. Hence 6 = 1/2. The solutions have been plotted at t = 1/2. 

u(x,0) - H(x,2,0), u{-1 ,t) - 2. 1-0.5. n-100 



Fig. 5: Moving Shock, 4th-order Explicit, 6—1/2 


u(x,0) « H(x,2.0), u{-1,t) - 2. 1-0.5. n-100 



Fig. 6: Moving Shock, 4th-order Implicit, £ = 1/2 

Next we employ the viscosity based on perturbation analysis as described in section 
3.2. The data is the same as in the previous examples. The viscosity is now given by 
eq. (47), i. e., t\ — 1 /3, £2 = —1/12. We also give two examples of the averaging technique 

of section 3.3. 
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We conclude this section by showing two 4th-order computations where we have solved 
the linear advection equation u t + u x = 0 by means of the hybrid scheme (61). In the first 
example we have used £ = 1/2 and k = h. The switch r, has been set to one at every 
grid point. The hybrid scheme is thus equivalent to the method of characteristics. In the 
second example £ = 1/2 and k = 0.9 fc; the switch rj is now turned on adaptively. This 
implies that the numerical method should behave as an approximate upwind scheme near 
the discontinuity. The solutions have been plotted at t = 1/2. 


ti 1 r 

8- 

6- 

4- 

2- 

r 1 * 

* 

, , , ... i — i 1 1 

-1 -0.8 -0.6 -0.4 -0.2 0 0.2 0.4 06 0.8 1 

a Discontinuity, 4th-order Explicit (+), — - — 

u(x,0) ■F(x,l), u(-1,t)-0. t • 0.5, n-200 

Z ' ■■ 


1 - 
8- 
6- 
4- 

.2- 

* 

1 i 

r 

-oft -0 6 -0 4 -0.2 0 0.2 0.4 0.6 0.6 1 


Fig. 11: Propagation of a Discontinuity, 4th-order Explicit (+), — • — Initial Data, k 
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